\documentclass[10pt]{article}
\usepackage[UTF8]{ctex}
\usepackage[top=2.5cm, bottom=2.5cm, left=2.5cm, right=2.5cm]{geometry} % 页边距
\usepackage{amsmath, amssymb, amsthm}
\usepackage{graphicx}
\usepackage{float}
\usepackage{hyperref}
\usepackage{listings}
\usepackage{xcolor}
\usepackage{verbatim}

\usepackage{titling}
\setlength{\droptitle}{-2cm} % 标题上移

\title{有空气阻力的单摆运动建模}
\author{AI ET AL}
\date{\today}

\begin{document}

\maketitle

\abstract{对有空气阻力的单摆运动进行常微分方程建模。1、设计实验，测量数据，    
使用模拟数据拟合参数；2、对常微分方程线性化，在角度-角速度(x-v)相图    
里分析原点（平衡点）的稳定性；3、根据所得模型进行仿真，画出单摆运动    
的动画。%将相关结果分别写在新的tex文件和py文件里。
}

\section{引言}

单摆是物理学中的经典模型，其运动方程在无阻尼情况下为:
\[
\frac{d^2\theta}{dt^2} + \frac{g}{l} \sin\theta = 0
\]

然而，在实际环境中，单摆会受到空气阻力的影响。本文将建立包含空气阻力的单摆模型，并进行参数拟合、稳定性分析和仿真。

\section{有空气阻力的单摆模型}

\subsection{物理模型建立}

考虑一个质量为$m$的小球，用长度为$l$的无质量刚性杆悬挂。设$\theta(t)$为摆角（从竖直向下为零，逆时针为正）。

作用在小球上的力包括：
\begin{enumerate}
  \item 重力：$mg$，切向分量为$-mg\sin\theta$
  \item 空气阻力：与速度方向相反，大小与速度成正比，为$-k v = -k l \frac{d\theta}{dt}$，其中$k$为阻尼系数
\end{enumerate}

根据牛顿第二定律，切向方向的运动方程为：
\[
ml\frac{d^2\theta}{dt^2} = -mg\sin\theta - k l \frac{d\theta}{dt}
\]

整理得：
\[
\frac{d^2\theta}{dt^2} + \frac{k}{m} \frac{d\theta}{dt} + \frac{g}{l} \sin\theta = 0
\]

这是一个二阶非线性常微分方程。

\subsection{参数说明}

\begin{itemize}
  \item $g = 9.8\, \text{m/s}^2$：重力加速度
  \item $l$：摆长（米）
  \item $m$：小球质量（千克）
  \item $k$：空气阻力系数（kg/s）
\end{itemize}

\section{模型线性化与平衡点分析}

\subsection{一阶系统转换}

令$x_1 = \theta$, $x_2 = \frac{d\theta}{dt}$，则原方程可转换为一阶系统：
\[
\begin{cases}
\frac{dx_1}{dt} = x_2 \\
\frac{dx_2}{dt} = -\frac{k}{m} x_2 - \frac{g}{l} \sin x_1
\end{cases}
\]

\subsection{平衡点}

令$\frac{dx_1}{dt} = 0$, $\frac{dx_2}{dt} = 0$，得平衡点：
\[
x_2 = 0,\quad \sin x_1 = 0 \Rightarrow x_1 = n\pi, \quad n \in \mathbb{Z}
\]

主要平衡点有两个：
\begin{enumerate}
  \item 稳定平衡点：$(0, 0)$（摆垂直向下）
  \item 不稳定平衡点：$(\pi, 0)$（摆垂直向上）
\end{enumerate}

\subsection{线性化分析}

在平衡点$(0, 0)$附近，由于$\sin x_1 \approx x_1$（小角度近似），系统线性化为：
\[
\begin{cases}
\frac{dx_1}{dt} = x_2 \\
\frac{dx_2}{dt} = -\frac{k}{m} x_2 - \frac{g}{l} x_1
\end{cases}
\]

写成矩阵形式：
\[
\frac{d}{dt}
\begin{pmatrix}
x_1 \\
x_2
\end{pmatrix}
=
\begin{pmatrix}
0 & 1 \\
-\frac{g}{l} & -\frac{k}{m}
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2
\end{pmatrix}
\]

设$A = \begin{pmatrix} 0 & 1 \\ -\frac{g}{l} & -\frac{k}{m} \end{pmatrix}$为线性化系统的雅可比矩阵。

\subsection{稳定性分析}

特征方程为：
\[
\det(A - \lambda I) = \lambda^2 + \frac{k}{m} \lambda + \frac{g}{l} = 0
\]

特征根为：
\[
\lambda = \frac{-\frac{k}{m} \pm \sqrt{(\frac{k}{m})^2 - 4\frac{g}{l}}}{2}
\]

记$\alpha = \frac{k}{m}$（阻尼比），$\omega_0^2 = \frac{g}{l}$（固有频率的平方）。

则特征方程为：$\lambda^2 + \alpha\lambda + \omega_0^2 = 0$

根据阻尼比$\alpha$与$2\omega_0$的关系，系统有三种响应模式：

\begin{enumerate}
  \item \textbf{过阻尼}：当$\alpha > 2\omega_0$时，两个不同的负实根，系统无振荡地回到平衡位置。
  \item \textbf{临界阻尼}：当$\alpha = 2\omega_0$时，两个相等的负实根，系统最快回到平衡位置而无振荡。
  \item \textbf{欠阻尼}：当$\alpha < 2\omega_0$时，一对具有负实部的共轭复根，系统振荡衰减至平衡位置。
\end{enumerate}

在所有情况下，由于特征根的实部均为负，平衡点$(0, 0)$都是渐近稳定的。

\subsection{相图分析}

在$x-v$平面（角度-角速度平面）中，系统的行为表现为：
\begin{itemize}
  \item 过阻尼：轨迹单调地趋向原点
  \item 临界阻尼：轨迹最快地趋向原点
  \item 欠阻尼：轨迹以螺旋形式趋向原点
\end{itemize}

这表明无论阻尼大小如何，阻尼力都会使系统最终趋向于稳定平衡点$(0, 0)$。

\section{参数拟合实验设计}

为了拟合参数$k$，我们设计如下实验：

\begin{enumerate}
  \item 使用相同质量$m=0.1\,\text{kg}$和摆长$l=1.0\,\text{m}$的小球。
  \item 给予小球初始角度$\theta_0 = 0.3\,\text{rad}$（约17度），初始角速度为0。
  \item 记录单摆角度随时间的变化数据。
  \item 使用模拟数据拟合阻尼系数$k$。
\end{enumerate}

\section{仿真与动画}

基于建立的数学模型和拟合的参数，我们可以进行数值仿真并生成单摆运动的动画。动画将展示：

\begin{enumerate}
  \item 单摆的实时运动轨迹
  \item 阻尼对摆动幅度的影响
  \item 系统趋向稳定平衡点的过程
\end{enumerate}

通过动画仿真，我们可以直观地观察到阻尼系统的行为特性，验证理论分析的正确性。

\section{代码说明}

\begin{verbatim}
   1. pendulum_config.py - 包含物理常数和微分方程定义
   2. generate_simulated_data.py - 生成模拟实验数据的函数     
   3. fit_model.py - 参数拟合函数
   4. plot_phase_diagram.py - 相图绘制函数
   5. create_animation.py - 动画生成函数
   6. main.py - 主程序，协调所有功能
\end{verbatim}

\begin{thebibliography}{9}
\bibitem{dingtongren} 丁同仁,李承治.常微分方程教程[M].第3版.高等教育出版社,2022.
\bibitem{sishoukui-2} 司守奎,孙玺菁.Python数学建模算法与应用[M].国防工业出版社,2023.
\end{thebibliography}

\end{document}